Accepted for publication in ApJ 

Preprint typeset using L^-T^X style cmulatcapj v. 5/2/11 



DIRECT NUMERICAL SIMULATION OF RADIATION PRESSURE-DRIVEN TURBULENCE AND WINDS IN 

STAR CLUSTERS AND GALACTIC DISKS 



Department of Astronomy 



Mark R. Krumholz 

Astrophysics, University of California, Santa Cruz, CA 95064 USA 



o 

(N 
O 

O 

o 

o 
u 

6 

CO 

(N 
> 

(N 

0\ 
(N 

en 

o 

(N 



* 



Todd A. Thompson 

Department of Astronomy and Center for Cosmology & Astro-Particle Physics 
The Ohio State University, Columbus, OH 43210-1173 USA 
Accepted for publication in ApJ 

ABSTRACT 

The pressure exerted by the radiation of young stars may be an important feedback mechanism that 
drives turbulence and winds in forming star clusters and the disks of starburst galaxies. However, 
there is great uncertainty in how efficiently radiation couples to matter in these high optical depth 
environments. In particular, it is unclear what levels of turbulence the radiation can produce, and 
whether the infrared radiation trapped by the dust opacity can give rise to heavily mass-loaded winds. 
In this paper we report a series of two-dimensional flux-limited diffusion radiation-hydrodynamics 
calculations performed with the code ORION in which we drive strong radiation fluxes through columns 
of dusty matter confined by gravity in order to answer these questions. We consider both systems 
where the radiation flux is sub-Eddington throughout the gas column, and those where it is super- 
Eddington at the midplane but sub-Eddington in the atmosphere. In the latter, we find that the 
radiation-matter interaction gives rise to radiation-driven Rayleigh- Taylor instability, which drives 
supersonic turbulence at a level sufficient to fully explain the turbulence seen in Galactic protocluster 
gas clouds, and to make a non-trivial contribution to the turbulence observed in starburst galaxy disks. 
However, the instability also produces a channel structure in which the radiation-matter interaction 
is reduced compared to time-steady analytic models because the radiation field is not fully trapped. 
For astrophysical parameters relevant to forming star clusters and starburst galaxies, we find that this 
effect reduces the net momentum deposition rate in the dusty gas by a factor of ~ 2 — 6 compared to 
simple analytic estimates, and that in steady state the Eddington ratio reaches unity and there are 
no strong winds. We provide an approximation formula, appropriate for implementation in analytic 
models and non-radiative simulations, for the force exerted by the infrared radiation field in this 
regime. 

Subject headings: galaxies: ISM — galaxies: star clusters — hydrodynamics — instabilities — ISM: 
jets and outflows — radiative transfer 



1. INTRODUCTION 



The idea that the pressure exerted by stellar ra- 
diation might have significant dynamical effec ts on 
the gas in galaxies is an old one (O'dcll ct al. 1967; 
Chiao fe Wickramasinghdll972tlElmegreenlll983l : lFerraral 



19931 : IScoville et all 120011 : IScovilldl^nni rbut it has re- 
ceived significant renewed attention in recent years as a 
possible explanation for various pheno mena in star clus- 
ters a nd g alaxies. On ga l actic scales, rThompso n et al.l 
(|2005ft and iMurrav et al.l (|2005l ) propose that the force 
exerted by radiation from newly-formed stars both sets 
an upper limit on the star formation rate in ultralumi- 
nous infrared galaxies (ULIRGs) and drives the highly 
supersonic motio ns that are ub iquitousl y obs erved in 
these systems. lAndrews fc Thompson! ([20111 ) extend 
this analysis and propose that radiation pressure lim- 
its on the star fo r matio n rat e even in normal ga lax- 
ies. IMurrav et all ()2005f ) and IMurrav et all (J2011I ) ar- 
gue that radiation pressure is also responsible for launch- 
ing galactic winds (se e also Zhang fc Thompson 2012). 
On subgalactic scales, [Krumholz & Matzncr (2009]) and 
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iFall et all ([20101 ) argue that radiation pressure is the 
dominant feedback mechanism for massive young star 
clusters, and that winds driven by radiation momentum 
set the star format ion efficiency i n clust ers and the clus- 
ter mass function. M urray et al.l (|2010f) also argue that 
radiation pressure is the primary feedback mechanism in 
massive star clusters, and that it is responsible for limit- 
ing the star formation efficiency of giant molecular clouds 
across a wide range of galactic environments at low- and 
high-redshift. 

This renewed theoretical attention has given rise to 
a number of approximate implementations of radiation 
pressure feedback in simulations of galaxy evolution. The 
e arliest of these is the "momen tum-driven wind" model 
of Oppcn heimer fc Pavel ([20061 ). in which they inject mo- 
mentum into star-forming gas in their cosmological sim- 
ulations at a rate that depends on both to the star for- 
mation rate and the depth of the galactic potential well; 
the latter dependence is an attempt to approximate the 
increase in radiative momentum imparted to gas that oc- 
curs when the optical depth is high, so every photon is 
absorbed and r eemitte d multip le times. More recently, 
iHopkins et"aH ([2011allb1 . I2012bllal ) have implemented a 



more sophisticated model for radiation feedback into 
their isolated galaxy simulations. In their approach, the 
code identifies contiguous star-forming clumps and then 
applies an outward radiation force to the gas in them. 
The force is proportional to the product of the luminosity 
produced by all stars in the clumps and the gas column 
density; as with the momentum-driven wind model, the 
latter dependence is an attempt to capture the effects of 
force amplification due to radiative trapping. 

However, both the analytic models and the resulting 
semi-analytic implementations of feedback in the simu- 
lations contain significant uncertainties. At a basic level, 
it is unclear whether the radiation is actually able to 
drive winds from gravitationally bound galaxies or gas 
clumps. The opacity of dusty material varies with tem- 
per ature_js_^cjaghlv_K ex T 2 (at temperatures < 150 
K; iSemenov et all 12003). and the high optical depths of 
dusty ULIRGs or star-forming clumps ensures that their 
central temperatures are higher than the temperatures 
at their edges. For such objects it is often the case that 
the Eddington ratio (defined as the ratio of the radiative 
and gravitational forces per unit mass) is larger than 
unity for the warm material near the center, but less 
than unity for the cooler material near the edge. In this 
case it is unclear whether ra di ation will laun c h a w ind at 
all. iThompson et all (120051) iMurrav et all ([2010D . and 
lAndrews fe Thompson! ( 20111 ) argue that such a config- 
uration will produce a wind if the temperature and thus 
the Eddington ratio at the cent er is sufficiently high, 
while Krumholz fc Matzncr (2009) argue that a wind will 
occur only if the Eddington ratio exceeds unity at the 
edge of the object. 

If there is a wind, a second uncertainty is in how much 
momentum the radiation deposits in the matter. Con- 
sider a source of radiation of luminosity L, such as a 
young star, surrounded by dusty gas. If every photon 
is absorbed once, the radiation will deposit its full mo- 
mentum flux L/c in the gas. However, if the medium 
is so optically thick that each photon is absorbed and 
re-emitted multiple times, the amount of momentum de- 
posited in the matter could be significantly larger. In 
the limit where every photon is absorbed many times, 
all the energy of the radiation field is transformed into 
kinetic energy of the gas, and the momentum transfer 
approaches L/v, where v is the gas characteristic veloc- 
ity, and the exact prefactor will depend on the gas ve- 
locity and density distribution. The two limiting cases 
of L/c and L/v may be referred to as the momentum- 
driven and energy-driven limits, respectively, since in the 
former case the momentum deposited is limited by the 
momentum of the radiation field, while in the latter case 
it is limited by the energy of the radiation field. 

Different authors have come to differing conclu- 
sions a bout where between t hese l imits systems will 
fall. Krumholz fc Matznerl (|2009j) argue that the 
momentu m deposited will n ever exceed a few L/c, 
while [T homps on et all (120051). IMurrav et all (|2010[ ). and 
lAndrews fc Thompson! (J2011I ) argue that in optically 
thick systems it will be ttrL/c, where tir is an appropri- 
ately defined mean IR optical depth, which can be ^> 10 

1 Note that L/c -C tjrL/c <C L/v, so models with tir are 
intermediate between the pure momentum- and energy-driven lim- 
its. However, in the literature feedback models with T[r are still 



Which argument turns out to be correct has important 
implications for questions like whether radiation pressure 
is the domi nant mechanism for disrupting most molec- 
ular clouds ([Murray et all 120101 ) or only those fo rming 
the most luminous star clusters ([Fall et alJ l2010l). and 
wheth er giant clumps in z ~ 2 galaxies survive for long 
times ([Krumholz fc Dekel 2 01~0l) or are rapi d ly disrupted 
by st ellar feedback ([Hopkins et all l2011at iGenel et all 
1201211 

The strength and nature of radiation-matter coupling 
is uncertain because the matter distribution in real galax- 
ies and star clusters is highly non-uniform, leading to 
complex density fields through which the radiation must 
pass. Moreover, there are numerous instabilities that 
can occur when radiation exerts strong forces on matter, 
such as the photon bubble instability ( Blaes fc Socrates! 
2003) and the radiation Raylcig h- Taylor insta bility 
([Krumholz et all 120091: Dacquet fc Krumholzll2011l ). and 
these alter the distribution of both matter and radia- 
tion. The question of radiation-matter coupling there- 
fore requires modeling the fully non-linear development 
of radiation- hydrodynamic instabilities, which in turn re- 
quires simulations capable of treating both the radia- 
tion field and the matter. Both calculations of radia- 
tive transfer through fixed density fields and calcula- 
tions of the density field in which the radiation field is 
taken to be un i form ( as is assumed, for example, in the 
Hopk ins et all (|2011b[ ) feedback prescription) are inad- 
equate to the task. To date no simulations capable of 
answering this question have been reported. 

In this paper we explore a simple model system that 
nonetheless contains many of the essential features re- 
quired to study non-linear matter-radiation coupling. 
We use this system to understand the nature of this cou- 
pling, and to derive estimates for current unknowns such 
as the ability of radiation pressure to drive turbulence 
and winds, and the efficiency with which radiation de- 
posits momentum in the gas. In Section [5] we describe 
our model system and the equations that govern it, ob- 
tain its important dimensionless numbers, and determine 
under what circumstances it has an equilibrium state. In 
Section [3] we describe the numerical methods we use to 
simulate our model system, and in Section 0] we describe 
the results of our numerical simulations, the limitations 
of our calculations, and caveats to our conclusions. Fi- 
nally, in Section [5] we discuss the implications of our re- 
sults, and in Section [6] we summarize. 

2. MODEL SYSTEM 

2.1. Description and Basic Equations 

We will treat a section of a galactic disk or a young 
star cluster as an idealized model system, which allows us 
to isolate the physics of the radiation-matter interaction 
without worry about the complexity of a real disk or 
cluster. We consider a slab of gas with total surface 
density £ filling the domain z > 0. It is subject to a 
constant gravitational force per unit mass — gz that is 
independent of the position x. A vertical radiation flux 
F = Foz enters the domain of interest at z = 0. We 
neglect the self-gravity of the gas and assume that all 

referred to as momentum-driven to distinguish them from mod- 
els in which the energy of hot supernova-shocked gas dominates 
feedback. 



radiation is injected at z = (i.e. there are no internal 
sources of radiation at z > except the thermal emission 
of the gas itself.) 

For simplicity, we adopt the two-temperature flux- 
limited diffusion approximation. This allows us to in- 
terpolate between the optically thin and optically thick 
limits approximately without the need to track the radi- 
ation spectrum at each point, while not forcing the gas 
and radiation temperatures to be reach equality instan- 
taneously in low optical depth regions where the matter 
and radiation are weakly coupled. In this approximation, 
the radiation flux F and energy density E are related by 
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where up is the Rosseland mean opacity (which can in 
general be a function of the gas temperature, the radi- 
ation energy density, and the gas density), p is the gas 
density, and A is the dimensionless flux limiter, given in 
detail below. The equation s of radiation hydrody namics 
applicable to this case are (Krumho lz et al.H2007T> 



m p = - V-(pv) 



(2) 



- (pv) - - V • (pvv) - VP - XV E - pgz (3) 

— (pe) = -V • [{pe + P)v] - K P p{^B - cE) 



X (2— -1 ) v- YE 
V K R 

pgv z 



3 - i? 2 v 2 



-Kpp — E 

2 c 



(4) 



— E = V- (—YE) + KppttnB - cE) 
at \ krp ) 

-\U^-\\vYE 

3-i?2 V 2 {3-R 2 

+ —-Ik p p-E-Y-{—-^vE). (5) 



where v is the gas velocity, T g is the gas temperature, 
P = phpTg/p is the gas pressure, e = (7 — l)~ 1 fcsX' g //j + 
v 2 /2 is the gas specific energy, p is the mean mass per gas 
particle, 7 is the gas ratio of specific heats, B = caT g /Air 
is the frequency-integrated Planck function, Kp is the 
Planck mean opacity, and R2 is the Eddington factor. 
For convenience we also define the radiation temperature 
by T r — (E/a) 1 ' . The above equations are written in 
the mixed frame, so Kp and Kp are to be evaluated in 
the frame comoving with the gas but all other quantities 
are evaluated in the lab frame, so that total energy is 
conserved. The functions A and R2 have the property 
that A — > 1/3 and R2 — > 1/3 in the optically thick limit, 
and A — > KppE/YE and R2 — > 1 in the optically thin 
limit. We adopt the flux limiter and Eddington factor of 
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A =l( COthi? -| 



R= 



\YE\ 



(G) 

(7) 
(8) 



K R pE 
R 2 = X + X 2 R. 

2.2. Non-Dimensionalization 

In order to characterize the problem better, and un- 
derstand what dimensionless numbers control it, we now 
non-dimensionalize equations © ^ ©. We first note 
that, in any steady state configuration, the energy enter- 
ing the slab at z = must match the energy escaping 
to z — 00. As z —¥ 00 the density must approach for 
any physically reasonable configuration, so for some suf- 
ficiently large z the gas must be optically thin, and thus 
as z — > 00 the flux and radiation energy density must 
approach the relationship F^ = cE^z. This motivates 
us to non-dimensionalize the equations by defining a ref- 
erence temperature 
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(9) 



In steady state, both T r and T g must approach T* as 
z — > 00. 
From this temperature we also define 
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as the associated sound speed, scale height, and density; 
here p, is the mean mass per particle in hydrogen masseso 
With these definitions, we make a change of variables 

£ = x//i* s = t/t* t* = h*/c s ^, (11) 

and equations dSJ — (EJ become 
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2 Note that this scale height is the value we would expect in the 
limit of negligible radiation forces, but even if ft,* does not describe 
the actual gas configuration, it still provides a convenient reference 
length. 
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where V^ indicates spatial differentiation with respect to 
£ rather than x, and we have defined the non-dimensional 
variables 
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(17) 
(18) 
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These are, respectively, the dimensionless density, ve- 
locity, gas temperature, radiation temperature, Planck 
mean opacity, Rosseland mean opacity, and gas specific 
energy. The dimensionless ratio R that determines the 
flux-limiter is given by 
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The dimensionless ratios appearing in equations (ji"2"j) - 
flUD are 
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K>R,* 

where kr^ = Kp(/9*,T») and similarly for ftp,*. These 
quantities have simple physical interpretations: for mat- 
ter and radiation at the reference density and temper- 
ature /3* and T«, /e,* is the ratio of the radiative and 
gravitational forces, p B is the ratio of the sound speed to 
the speed of light, r* is the Rosseland mean optical depth 
of the slab of gas, and fco is the ratio of the Planck and 
Rosseland mean opacities. Finally, the condition that 
the flux at z = be Fq combined with the flux-limited 
diffusion approximation (equation [T]) requires that 



dQ r 
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at z = 0. 

An additional simplification is possible if we limit our- 
selves to the static diffusion or streaming limits, meaning 



that we require that /3 s t* <C 10 This is likely to hold 
in any real galactic disk or star cluster, since r* is never 
more than a few tens, while j3 s is generally of order 10~ D . 
In this case we can drop the terms proportional to /3 s r* 
in equations (|14|) and (|15[) . and these equations simplify 
to 

y s (b q ) = -v ? [b ( q + e g ) u] - ^h^k k P b(e 4 g - e*) 

+ fel\ ( 2k ^- - 1 ) u • Ve& 4 r - bu 7 (26) 
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The above analysis demonstrates that for a given func- 
tional form of the dimensionless opacities kp and fcp, the 
behavior of the fluid in this simplified set up is dictated 
by only three dimensionless numberfl the optical depth 
r*, the Eddington ratio /e,*, and the dimensionless ratio 
(3 S . The first of these determines how effectively radiation 
is trapped in the slab, the second determines the dynam- 
ical importance of radiation pressure relative to gravity, 
and the third determines the relative importance of ra- 
diation advection to radiation emission and absorption 
in matter radiation coupling, though its precise value is 
unlikely to matter as long as f3 s <C 1. Thus in practice 
t* and /e,* determine the parameter space of interest. 
Given values for the dimensionless parameters, the di- 
mensionless solution can then be scaled to physical units 
by a choice of the mean particle mass p and the two 
dimensional parameters £ and g. 

2.3. Equilibrium Solutions 

To understand the behavior of the equations, it is help- 
ful to first search for equilibrium solutions, in which all 
time derivatives vanish and u = 0. Examination of equa- 
tion (|14p immediately indicates that such solutions must 



obey O r = S = ©. Inserting this condition and u = 
into the remaining equations reduces the problem to a 
pair of coupled, nonlinear ordinary differential equations 



^_(6e) + 4A Ae 3^ + 6 = o 
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(28) 
(29) 



Note that this system of equations depends only on /e,* 
and r*, not on j3 s or fco- The latter two quantities are 
relevant only when the radiation field is out of equilib- 
rium. These equations are to be integrated from z = 

3 Formally, the product /3 s t„ determines whether the system is 
described by static or dynamic diffusion. For more discussion, see 



IKrumholz et ail f2007T) . 

4 The opacity ratio fco also enters, but it is always of order unity 
for physically reasonable continuum opacity sources. 



to oo, subject to the boundary conditions 
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which are equivalent to requiring that the flux be Fq 
at z — and z — oo. The third boundary condition 
required to specify this third-order system conies from 
the integral constraint that 



6d£, = 1, 



(32) 
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which is equivalent to demanding that J pdz — S. 

We can analytically integrate equation (|2T)]) once, and 
doing so and using the boundary condition (|30p allows us 
to rewrite the system in the somewhat more transparent 
form 



d& 



These equations have a simple physical interpretation. 
The quantity 60 is the dimensionless gas pressure, and 
/e,*&r is the dimensionless version of the Eddington ra- 
tio krF '/ ' gc at a given point in the disk. Since the flux is 
invariant with z, this in turn is just the Eddington ratio 
at infinity, /e,*, scaled by kn, the ratio of the local opac- 
ity to the opacity at infinity. Thus equation (1331) simply 
asserts that the gas pressure gradient balances the force 
of gravity, diluted by radiation pressure, at every point. 
Equation (|34[) asserts that the temperature gradient is 
such that the radiation flux is constant with height. 

For a given value of /e,* and t„, and a specified func- 
tional form of fcfl, we can solve equations (|33p and (|34l) 
via a double iteration procedure. We begin by guessing 
values for b(0) (which must be < 1) and 0(0) (which 
must be > 1) at z = 0, and using these boundary con- 
ditions to integrate the system of ODEs from z = up 
to a value of z large enough so that d®/dt; z w 0, or until 
we encounter a singular point where — > at finite z. 
In general the resulting solution will not satisfy the con- 
straint that — > 1 at large z (equation [3"Tj) . We therefore 
hold b(Q) fixed and iterate on 0(0) until we find the value 
of 0(0) that does give — > 1 at large z. However, this 
will still not generally satisfy the integral constraint on 
b (equation I32[) , and we must therefore iterate again to 
find a value of 6(0) such that equation (|3"2")l is satisfied. 
For each value of 6(0) we must again iterate on 0(0), 
giving rise to a double iteration^ 

Figure [1] shows two sample solutions, both computed 
with kjf = 2 , the functional form that approximately 
descr ibes dust opacity in the infrared ()Semenov et al.l 
2003). The two solutions shown correspond to /e,* = 0.3, 
r* = 1 and /e,* = 0.03, t* = 10. These choices do not 
map to a unique set of physical parameters, but as an 

5 In principal we could eliminate one iteration by adopting = 1 
at some large z as a boundary condition and integrating in the — z 
direction. In this direction, however, the system of ODEs is stiff, 
and thus it is more computationally efficient to integrate in the +z 
direction and iterate on both b(0) and 0(0). 




Fig. 1. — Equilibrium values for b (blue) and (red) as a func- 
tion of £ z , computed for an opacity law kn = 2 , with /e,* = 0.3 
and T* = 1 (solid) and /e,» = 0.03 and t* = 10 (dashed). The 
Eddington ratios at § z = for these solutions are 0.82 and 0.49, 
respectively. 

example choice we can consider the same flux value and 
opacity law that we will use in our numerical experiments 
below. The flux is F = 2.6 x 10 13 L Q kpc~ 2 , appro- 
priate for a bright ULIRG; the corresponding value of 
T; = 82 K, and this gives k r .* = 2.1 cm 2 g" 1 . With 
this opacity and radiation flux, the remaining dimen- 
sional parameters are £ = (0.47,4.7) g cm~ 2 , g — (2.5 x 
10" 6 , 2.5 x 10~ 5 ) dyne g~\ h* = (3.8 x 10~ 4 , 3.8 x 10~ 5 ) 
pc, and p* = (4.0 x 10~ 16 , 4.0 x 10 -14 g cm -3 ), where the 
first numerical value given in parentheses is for r* = 1, 
/e,* = 0.3, and the second is for t* = 10, /e,* = 0.03. 

The qualitative behavior of the solutions is straight- 
forward to understand. At large £ z where the gas is op- 
tically thin, the equilibrium dimensionless temperature 
approaches a constant value = 1. However, the tem- 
perature at the midplane (£ z = 0) is higher because the 
gas is optically thick. The increase in temperature at the 
midplane is largest for the case t* = 10, since a higher 
optical depth leads to more effective trapping of the ra- 
diation. In both cases, the increased temperature leads 
to a decreased density 6 at £ z = compared to b = 1 , the 
value that would be produced if the gas were isothermal. 
In addition, the density declines much more slowly with 
£ z than would be the case for a simple isothermal atmo- 
sphere; the scale height is ~ 5 — 10, compared to 1 for a 
simple isothermal atmosphere. The more gradual falloff 
in the density arises from two effects. First, since the 
Eddington ratio near the midplane is a significant frac- 
tion of unity, radiation pressure force helps support the 
atmosphere against gravity. Second, the gas near the 
midplane is hotter than for a simple isothermal atmo- 
sphere, further increasing the scale height. Both of these 
effects are larger in the /e,* = 0.03, t* = 10 case due to 
the greater optical depth of the atmosphere, even though 
the radiation pressure force at the top of the atmosphere 
in this case is smaller than in the /e,* = 0.3, t* = 1 case. 
Both of these effects are strongest at small values of £ z 
where the temperature, radiation energy density, and ra- 
diation force are elevated. At larger £ z , where © sw 1, 
radiation pressure force becomes weak, the gas is close 
to isothermal, and the density falls off rapidly, returning 
to the behavior expected for a normal isothermal atmo- 
sphere. 

The iteration procedure we use to generate these so- 
lutions is not guaranteed to converge, because an equi- 




3.2. Choice of Simulation Parameters 
In all simulations we adopt an opacity law 

K(R , P) = (10- 3 /M0- 1 )^) 2 cm 2 g- 1 , 



(35) 



Fig. 2. — /E.criti the maximum value of /e,, at the slab surface 
for which equilibrium is possible, as a function of slab optical depth 
r«, using an opacity law kji = © 2 . Circles indicate values of r* for 
which we numerically determined values of /e cr ; t . 

librium profile for which the density and temperature 
remain finite as z — > oo is not guaranteed to exist for an 
arbitrary combination of /e,*, t*, and k^. In particular, 
note that if /e,* > 1, then since kn — > 1 as £ z — y oo, it 
follows that the right hand side of equation (|3"3")l is non- 
negative as £ 2 — » oo. In this case Jbd^ z diverges, and 
there is no solution possible that satisfies condition (|32j) . 
Moreover, even if /e.* < 1, there may still be no solution 
that obeys both the constraint equations (J3TJ) and (|32j) . 
In practice the maximum value of /e,* for which an equi- 
librium solution exists, which we refer to as /E,criti must 
be determined numerically for a given functional form of 
kfj and optical depth r* . Figure [2] shows the results of 
such a calculation for kp = O 2 . 

3. NUMERICAL SIMULATIONS 

Now that we have identified the important dimension- 
less numbers and found equilibria whenever they exist, 
we turn to a full numerical simulation that will allow us 
to explore the behavior of non-equilibrium systems. 

3.1. Numerical Method 

Our simulations use the OR ION radi ation- 
hydrodynamics code (|Kleinl 119991: iFisherl 1200% 
iKrumholz et al.l 12007ft : the properties of this code 
may be found elsewhere in the literature, so we simply 
summarize them here. ORION solves the equations of 
radiation-hydrodynamics (equations [5] - EJ) in the two- 
temperature flux-limited diffusion approximation. The 
code uses the conserv ative operator-splitting scheme of 
IKrumholz et al.l ()2007l ) to separate the implicit radiation 
and explicit hydrodynamic updates. The latter uses a 
high-order Godunov method that requires very little 
artificial viscosity (see |Kleinl 119991 fo r detai ls). The 
former uses the iShestakov fc Offnerl (|2008| ) pseudo- 
transient continuation method to solve the implicit 
radiation system. For the purposes of the simulations 
here, we do not use ORION 's additional capabilities for 
self-gravity and sink particles, and instead we impose a 
constant gravitational acceleration in the —z direction, 
per equation ([3]). Although ORION has adaptive mesh 
refinement capability, we do not use the AMR in these 
simulations, because we find that in most simulations 
the dense gas occupies a large fraction of the simulation 
volume for much of the computation time. This negates 
the computational advantage from using AMR. 



which is rou g hly i n accord with the model of 
iSemenov et al.l (|2003t ) at temperatures < 150 K. We 
choose not to use the full ISemenov et al.l opacity func- 
tion (which ORION includes) in order to keep the prob- 
lem as pure and simple as possible^ We note that this 
approximation will, if anything, lead us to overestimate 
the strength of matter-radiation coupling, since we are 
not including the flattening of the opacity at high tem- 
peratures. 

For our choice of other parameters, we pick values 
that span an interesting physical range, and that over- 
lap with observations. On the latter point, observa- 
tio ns of ULIRGs and the models to fit them provided 
bv lThompson et al.1 (|2005h give typical fluxes 10 13 - 10 14 
Lq kpc -2 , typical surface temperatures of T ~ 50 — 100 
K, and typical Rosseland mean optical depths at this 
temperature are t* ~ 1 — 10. The corresponding surface 
densities are ~ 1 — 10 g cm~ 2 , and the corresponding 
gravitational acceleration is g = 2ttGH ~ 10~ 6 — 10~ 5 
dyne g _1 , where we have computed g as for an infinite 
slab, and we have not included stellar mass. Combining 
these estimates we find values of /e,* in the range 0.01 — 1. 
If we focus on young clusters we find similar surface den- 
sities and fluxes, but with a som ewhat larger range, so 
that /e,* > 1 in some cases (e.g. IKrumholz fc Matznerl 
l2009f h Physically, the most interesting range to explore 
is the one where /E.crit < /e,* < 1. (Note that this 
regime exists only if r* > 1, since /E,crit is significantly 
less than unity only if the gas is optically thick) . These 
are the cases where there is no stable equilibrium, but 
the radiation force is not so strong that all the gas is 
certain to be ejected; instead the outcome will depend 
on how the radiation and matter couple. The majority 
of ULIRGs and some young, bright clusters are in this 
regime. 

Based on these considerations, we select a set of run 
parameters described in Table [1] These parameters have 
a range of optical depths r* and Eddington factors /e,*- 
We include one case in the stable regime, /e,* < /E,crit, 
as a control, and the remaining runs are in the unsta- 
ble, /e,* > /E,crit, regime. We do not include runs with 
/e,* > 1. Partly this is because the outcome in this case 
seems almost certain to be large-scale ejection, and partly 
because this regime does not appear to be reached in ob- 
served ULIRGs, although the latter statement is subject 
to significant uncertainties on the dust opacity and the 
CO to H2 conversion factor in ULIRGs, both of which 
affect estimates of /e,*. 

3.3. Simulation Setup 

All our simulations are two-dimensional Cartesian, 
taking place in the (x, z) plane. The force of gravity 

6 For numerical reasons we cap the opacity at the value corre- 
sponding to T = 10T*, and we set Kt R ^p\ = in material with 
density < 10~ 10 p*. These choices allow us to introduce a hot, 
zero-opacity ambient medium that can act as a constant pressure 
boundary condition at high altitudes. 



TABLE 1 
Simulation Physical Parameters 



Run Name 



Dimcnsionlcss Parameters 

T * /e,* /E,»//E,crit 



Dimensional Parameters 
£ g h* U 

[g cm" 2 ] [1CT 6 dyne g" 1 ] [10~ 4 pc] [kyr] 



T10F0.02 


10 


0.02 


0.51 


4.7 


37 


0.25 


0.045 


T03F0.50 


3 


0.5 


.3.8 


2.9 


3.7 


2.5 


0.46 


T10F0.25 


10 


0.25 


6.4 


4.7 


2.9 


3.2 


0.57 


T10F0.50 


10 


0.5 


12.8 


4.7 


1.5 


6.3 


1.1 


T10F0.90 


10 


0.9 


23.1 


4.7 


0.83 


11 


2.1 



Note. — All simulations use f3 s — 1.8 x 10 _tl and &q — 10 ' , which correspond to dimensional parameters 
T* — 82 K and Fa — 2.6 X 10 Lq kpc - . Dimcnsionlcss and dimensional parameters are related by using the 
opacity law given by equation J35D , and adopting a mean particle mass fi — 2.33, as expected for a gas of molecular 
hydrogen and helium mixed in the standard cosmic abundance. Values of /E,crit have been computed numerically; 
for t* = 3, / B( crit = 0.13, and for r* = 10, / E ,crit = 0.039. 



TABLE 2 

Simulation Numerical Parameters 



Run Name 


Ax/h* 


[L x x L„]/h* 


[N x x N z ] 


trun/t* 


T10F0.02 


0.5 


512 x 256 


1024 x 512 


78 


T03F0.50 


0.5 


512 x 1024 


1024 X 2048 


217 


T10F0.25 


0.5 


512 x 2048 


1024 x 4096 


210 


T10F0.50 


0.5 


512 x 2048 


1024 x 4096 


215 


T10F0.90 


0.5 


512 x 4096 


1024 x 8192 


205 


T10F0.25.LR a 


1.0 


512 x 2048 


512 x 2048 


276 



Note. — Col. 2: Grid spacing. Col. 3: Size of the computational 
domain in the x and z directions. Col. 4: Number of computational 
cells in the x and z directions. Col. 5: Time for which the simulation 
was run. 

a Physical parameters for this run arc identical to those given for run 
T10F0.25 in Table[T] 

is in the — z direction. The numerical resolution and 
other run parameters are described in Table [21 and we 
describe a resolution study we conducted to ensure that 
the resolution is adequate in the Appendix. Our bound- 
ary conditions are periodic in the x direction. At z = 0, 
we use reflecting boundary conditions on the hydrody- 
namic variables, and Neumann boundary conditions on 
the radiation field, with the radiation flux into the box 
set to the value corresponding to the parameters indi- 
cated in Table [1] At the upper z boundary, for our 
hydrodynamic boundary condition we place immediately 
outside the computational domain a low density ambient 
medium with density 10 -13 p*, velocity 0, and tempera- 
ture 10 3 T,|,; this places it in pressure balance in the ini- 
tial condition (see below), but allows dense matter with 
an upward velocity to leave the computational domain 
freely. For the radiation, we impose a Dirichlct boundary 
condition that the radiation energy density is E = aT^. 
We initialize all runs with a gas density distribution 

p = [1 + A p sin(27ra;/Ap)] 

(36) 
a temperature T — T*, and a velocity of 0, where A p and 
A p are the amplitude and wavelength of a perturbation 
we introduce to seed instability. In all simulations we use 
A p = 0.25 and X p = 256/i*. We initialize the radiation 
energy density to E = aT^ at all points. Thus in the 
absence of a radiation flux entering the domain at z = 0, 
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Fig. 3. — Time series showing the density (top) and radiation 
temperature (bottom) as a function of time in run T10F0.02. Gas 
temperatures are nearly identical to radiation temperatures every- 
where except where p/p* < 10~ 10 , where the set the opacity to 
zero. Note that the simulation domain extends to 256fe» in the 
vertical direction, but we show only the region from to 64h* . 

and for A p — 0, the region with density > 10 _10 p* would 
simply be a stable, isothermal atmosphere. Setting A p to 
a non-zero value but leaving the temperature and radia- 
tion field fixed would result in a series of stable gravity 
waves propagating through the atmosphere. 

4. SIMULATION RESULTS 

4.1. The Stable Case 

We first examine run T10F0.02 (recall our naming 
convention that T10 means r* = 10, and F0.02 means 
/e,* = 0.02), which we may think of as a sort of con- 
trol, in the sense that /E,*//E,crit < 1- This simulation 
therefore has a stable configuration toward which it is 
able to converge. Figure [3] shows a series of snapshots 
of the density and temperature distribution in the run. 
As the plot shows, the gas is initially pushed upward 
by radiation pressure. This is not surprising, because 
the initial configuration is an equilibrium appropriate to 
the case of an isothermal gas with no radiation pressure. 




Fig. 4. — Comparison between run T10F0.02 and the analytic 
solution for T* = 10 , /e .« = 0.02, computed via the procedure 
described in Section 12.31 The top panel shows density and the 
bottom shows radiation temperature; gas temperature is nearly 
identical. In both panels, the solid line is the analytic equilibrium 
solution, and asterisks show the mean density/temperature at a 
given vertical position in the simulation, with the mean computed 
over all horizontal positions and all simulations times t/t t > 70. 
The gray band shows the range of values found at a given vertical 
position over all times and horizontal positions. The dashed lines 
show the initial conditions in the simulation. 

Once radiation is injected into the simulation domain, 
the temperature near z = reaches ~ 4 — 5T*. Since the 
opacity varies as kr oc T 2 , this means that the opacity 
and the radiative force felt by this gas is increased by a 
factor of ~ 20 compared to gas at a temperature T* , and 
the local Eddington ratio in this gas approaches ~ 0.5. 
This accounts for the upward movement. Since radiation 
is trapped more effectively at horizontal locations where 
the initial density is somewhat higher, due to the initial 
perturbation, the upward motion is strongest there. This 
amplifies the initial perturbation, so that the horizontal 
variation visible at t/t* — 10 and 20 is much larger than 
that present in the initial conditions. After the initial 
transient however, the gas settles back down, and the 
horizontal fluctuations damp out. There is little change 
past t/t* ~ 50. 

In Figure 2] we compare the state of the system at late 
times to the equilibrium configuration for t» — 10 and 
/e,* = 0.02. As the plot shows, the density and tem- 
perature converge to the equilibrium solution very well. 
There is some oscillation at high altitudes as the system 
rings down, but the convergence is clear. It is worth not- 
ing that, even though /E,*//E,crit < 1 and / E ,» = 0.02, 
radiation pressure is non-negligible in the equilibrium 
configuration. Since kr oc T 2 , the peak temperature 
of T w 4T* at the base of the simulation domain corre- 
sponds to a local Eddington ratio /e ~ 1/3. 

4.2. Unstable Cases: Morphology and Overall Evolution 

We now turn to runs T10F0.25, T03F0.50, T10F0.50, 
and T10F0.90, for which no equilibrium solution exists. 
Figures [5] and H] show the time evolution of the gas and 
radiation energy density distributions in each of them. 
As in run T10F0.02, the gas is initially flung upward by 
radiation pressure. The gas is driven into a thin shell, 



with high radiation temperature beneath it and low ra- 
diation temperature above it. As time passes, however, 
the shell begins to buckle, developing fingers of gas that 
penetrate down into the radiation-dominated lower re- 
gion. The buckling is evident by £/£» ~ 5, fingers appear 
by t/t* ~ 10 — 20, and by t/t* ~ 50 a clear non-linear 
instability has set in. The dense fingers of gas reaching 
into the radiation-dominated region have much smaller 
upward velocities than the rest of the gas, and at late 
times they begin to fall back toward z — 0. 

As the instability develops the gas becomes turbulent. 
It also develops a clear channel morphology, with most 
of the gas mass in dense, nearly vertical filaments, and 
most of the volume filled by low-density gas between 
the filaments. Gas in the dense filaments is predomi- 
nantly falling rather than rising. As a result, in runs 
T03F0.50 and T10F0.25 none of the mass reaches the 
upper boundary of the computational domain. In runs 
T10F0.50 and T10F0.90 some mass does reach the top 
of the computational domain, but it is a small fraction 
of the total, and most of it is simply coasting from the 
initial launch at early times, rather than being actively 
driven upward. By t/t* ~ 100 — 150, the upward mo- 
tion has completely halted and the gas has fallen back 
to the midplane. Thereafter the system exhibits contin- 
uous turbulent motion with a roughly constant velocity 
dispersion and gas scale height. This appears to be the 
final, statistical steady state. 

4.3. Turbulence 

One of the important questions regarding radiation in- 
stabilities is whether they can explain, or at least con- 
tribute to, the large turbulent velocity dispersions seen 
in radiation-dominated galactic disks. To address this 
question we compute the mass-weighted horizontal and 
vertical velocity dispersions o~ x and a z as a function of 
time in our simulations. We define these by 
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(37) 
(38) 



where M is the amount of mass in the simulation domain 
and the integrals are over the full simulation volume. 

Figure [7] shows the result, and Table [3] summarizes 
what is shown in the Figure. Not surprisingly, in the 
stable run T10F0.02, after the initial transient the ve- 
locity dispersion is highly subsonic. It is dominated by 
the horizontal component, and oscillates up and down as 
the system rings down toward equilibrium. In the other 
runs we see that, after the initial transient, the veloc- 
ity dispersion approaches a roughly constant, supersonic 
value. The vertical component of the velocity dispersion 
is somewhat larger than the horizontal one. The highest 
Mach numbers we reach are ~ 10. 

4.4. Eddington Ratio and Trapping Factor 

Another quantity of interest is the net force applied by 
the radiation to the gas, and how this compares to the 
force of gravity. A closely related question is how much 
momentum the gas is able to extract from the radiation 
field and transfer into a wind. Recall that the radiation 
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Fig. 5. — Gas density distribution as a function of time in simulations T03F0.25, T10F0.25, T10F0.50, and T10F0.90, as indicated. Each 
row represents a simulation, and each column is at a fixed time, as indicated. Blank panels indicate that the simulation was halted before 
the indicated time. For run T03F0.25, the top of the computational domain only extends to 1024h», indicated by the dashed line. For run 
T10F0.90, the computational domain extends to 4096ft* , i.e. twice the vertical extent shown. Note that the dynamic range of densities in 
the simulation is much larger than the range 10~ 4 — 10 — *p* that we show. We pick this range because ~ 95% of the mass in the simulation 
volume lies within it at most times. 
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Fig. 6. — Same as Figure [5] but showing the radiation temperature r = T r /T*. The gas temperature distribution is very similar 
everywhere except at very high altitudes and low densities. 
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Run Name 



TABLE 3 

Simulation Outcomes 



CTx/Ca,* (Tz/Cs,* CTJCs.t (/e) /trap /tra 



K(T mp )S m Tin 



T10F0.02 


0.22 


0.13 


0.27 


0.18 


88 


35 


160 


T03F0.50 


3.2 


2.6 


4.1 


1.0 


5.0 


2.5 


15 


T10F0.25 


1.8 


3.3 


5.8 


1.0 


39 


25 


93 


T10F0.50 


6.5 


5.9 


8.8 


1.1 


22 


13 


120 


T10F0.90 


5.7 


17.1 


18.1 


1.2 


12 


7.0 


64 


T10F0.25.LR a 


1.2 


2.9 


5.1 


1.0 


10 


21 


110 



NOTE. — All quantities shown represent time averages over all times t/t. M > 175 in all runs 
except T10F0.02, where the average is over t/t* > 50. The quantity /trap is the trapping 
factor considering all material (as defined by cquation |42t , while /trap,w is the trapping factor 
computed considering only material with v z > 0. The quantity K(T mp )S RS tir. is the average 
optical depth at the end of the calculation, computed using the mass- weighted mean midplane 
temperature. (Volume-weighting gives a nearly identical result.) Values of tir at the start of 
the calculation are very similar. 
a Physical parameters for this run are identical to those given for run T10F0.25 in Table fll 





FlG. 7. — Velocity dispersion a versus time in runs T10F0.02, 
T03F0.50, T10F0.25, T10F0.50, and T10F0.90, as indicated in each 
panel. We show the x and z components u x and o z (thin dashed 
and solid lines, respectively), and the total velocity dispersion a = 
(al + v 2 z ) 1/2 (thick solid line). 



force per unit volume exerted on the matter is 
krpF 



fr 



:ad 



-WE, 



(39) 



where F is the radiation flux, and the second equality fol- 
lows from the flux-limited diffusion approximation. The 
mean vertical radiation force per unit area applied to the 
gas in the simulation domain is 



(/. 



rad,j 



l 



f rad • z dV, 



(40) 



where L x is the horizontal length of the computational 
domain, and the integral is over the full simulation 
volume. It useful to compare this force to two other 
quantities. One is the mean force exerted by gravity, 
f g = — T,g z . We define the mean Eddington ratio of the 



Fig. 8. — Mean Eddington ratio (/e) (left axis) and trapping 
factor /trap (right axis) as a function of time for runs T10F0.02, 
T03F0.50, T10F0.25, T10F0.50, and T10F0.90, as indicated in each 
panel. The solid line is the radiation force applied to all the mat- 
ter in the computational domain, and the dashed line is the force 
applied only to matter with v z > 0. The dotted horizontal line 
marks (/e) = 1. Note that all simulations begin with (/e) = 
because we initialize the simulation with no gradient in the radia- 
tion energy density; the rapid rise of (/e) to its initial plateau is 
the result of the radiation field reaching equilibrium rapidly, on a 
timescale that is not resolved in the plot. 

computational volume by 



</e> = 



\/rad,2/ 

Zg z 



(41) 



The bulk of the matter can be ejected by radiation pres- 
sure only if (/e) > 1, although a wind containing a small 
fraction of the matter may be driven even if (/e) < 1. 

The second useful comparison is between the force ex- 
erted by the radiation and the momentum flux carried 
by the radiation field, which is Fq/c. As discussed in 
Section [T] if the medium is optically thick enough to 
ensure that every photon is absorbed at least once (as 



12 



is the case for our simulations), we expect the radia- 
tion field to transfer at least this much momentum to 
the gas. However, the amount of momentum extracted 
could be significantly larger if every photon is absorbed 
many times, up to a maximum value of order Fq/v, where 
v is the gas characteristic velocity. To quantify where 
between these two limits ou r simulations fall, we follow 
iKrumholz fc Matznerl (|2009f ) and define the trapping fac- 
tor /t ra p by 



\/rad,z) — (1 + /trap) 



F 



(42) 



In the limit of one absorption per photon we expect 
/trap = 0, and for infinitely many absorptions we expect 
/trap ~ c/v. Models that assume strong trapping gen- 
erally adopt /trap ~ 7ia, where ttr is a fiducial optical 
depth that depends on the mean gas te mperature. For 
exampl e, in their numerical simulations IHopkins et al.l 
(|2011b[ ) adopt /t rap = max(0, kqY, — 1), where /to = 5 
cm 2 g _1 and £ is the gas surface density; values of KoT, 
in the simulations range from <C 1 up to ~ 50. Finally, 
note that (/e) and /t rap are related by 



(/ E )=(l + /trap)^. 



(43) 



We show (/ E ) and / trap for our simulations in Figure[51 
and summarize the mean values in Table [H The are sev- 
eral interesting points to be taken from this plot. First, 
in the stable run T10F0.02, (/ E ) < 1 at all times. This is 
not surprising, since the run parameters were chosen to 
have this feature. In contrast, in all of the unstable runs 
the pattern is that (/e) > 1 at early times before the in- 
stability sets in. Once the instability is established, (/e) 
drops to less than unity, and the gas falls back. Finally, 
at late times when the instability is in steady state, (/e) 
is very close to unity. The value of (/e) in this steady 
state does not depend on /e.* or t». Thus we see that, 
even though t*/e,* > 1 in all three of the unstable runs, 
the actual force applied to the matter by the radiation 
field saturates at the Eddington limit. Thus there is no 
large scale ejection of matter, despite the fact that /Edd 
is much larger than unity at the midplane at the start of 
the calculation. 

A similar effect is apparent if we consider the trapping 
factor /trap- In the stable run, T10F0.02, / trap - 100. 
This is because the radiation field is effectively trapped 
by the uniform matter distribution, so every photon is 
absorbed and re-emitted many times. In contrast, in 
the unstable runs /trap is significantly smaller once the 
instability develops. Roughly half this force is applied to 
material with an upward velocity, simply because roughly 
half the gas has a positive velocity at any given time. 
However, the identity of individual upward-moving and 
downward-moving pockets of gas is continually changing, 
so no material has a sustained positive velocity and is 
launched into a wind. A reasonable conjecture based on 
the simulation results is that in the regime where /e.* < 
1 < t */e,*, the gas self-adjusts to ensure (/e) « 1, and 
that /trap goes to whatever value is required to make 
this happen. In this case we would predict that /trap ~ 
t*//e,* — 1- However, all of this momentum is delivered 
to bound gas, and does not produce a wind. 



4.5. Conversion of Radiative Power to Kinetic, 
Thermal, and Gravitational Potential Energy 

A final quantity of interest is how much radiative power 
is converted to the kinetic, thermal, and gravitational 
potential energy of the gas. At late times, once the tur- 
bulence has reached steady state, the net power transfer 
rate between the various energy reservoirs must be zero. 
This conclusion follows simply from the fact that a, the 
mean gas temperature, and the gas scale height all ap- 
proach statistically constant values at late times, imply- 
ing that the kinetic, thermal, and gravitational potential 
energies of the gas must be statistically constant as well. 
In practice the system produces this effect via a balance 
between two processes. The radiation exerts forces on the 
gas, doing work; this effect is represented by the equal 
and opposite terms proportional to v and E in Equations 
@ and ([5]). Once the gas becomes turbulent, however, 
regions of compression appear, and in these regions ki- 
netic energy is converted to thermal energy, which in turn 
is turned back into radiative energy, a process described 
by the terms ±K,pp(AnB — cE) in Equations (j4} and §5§ . 
(Both exchanges can go in the opposite direction too: 
gas can do work on the radiation field, and cool gas can 
be radiatively heated. In our problem these effects are 
generally smaller than their opposites, however.) These 
rates of exchange are on average equal and opposite, and 
the level of turbulence self-adjusts to ensure that they 
remain so. 

In this respect a system that does not drive a wind, as 
we find in all our simulations, is fundamentally different 
than one that does drive a wind. In a system with a wind, 
the baryons that reach infinity carry kinetic, thermal, 
and gravitational potential energy, and thus some non- 
zero fraction of the radiant energy must be converted to 
other forms in the process of driving the wind. Without 
a wind, once the system reaches steady state the time- 
averaged net conversion of radiant energy to other forms 
is necessarily zeroQ This analysis leaves open, however, 
the related question of what level of turbulence must be 
created so that this steady state is achieved. Clearly 
the equilibrium value of a/c Si * must be a function of r, 
and /e,*- In this work we have sampled the (t»,/e,*) 
plane only sparsely, and so we are not yet in a position 
to construct a model for this mapping. However, we plan 
to revisit this topic in future work. 

5. DISCUSSION 

5.1. What is the Nature of the Instability? 

The instability that appears in our simulations 
is almost certainly the radiation Rayleigh- Taylor 
(RRT) instabili ty , firs t described in simulations by 
IKrumholz et al.1 (|2009f ) and formally analyzed by 

7 One might object that in a planar geometry such as ours a 
wind can never occur, since the gravitational potential energy of a 
gas parcel diverges at large £ z , unlike in spherical geometry. This 
is true in principle bu t turn s out to be irrelevant in practice. As 
we discuss in Section 15.5.21 inserting realistic astrophysical scal- 
ings into our dimensionless problem shows that the gas in our 
simulations reaches such small heights that curvature effects are 
negligible, and a planar approximation is in fact valid. Gas in our 
simulations is not launched into a wind not because we have placed 
it in a gravitational well that is infinitely deep, but because the ra- 
diation force is weak enough that it could not escape a potential 
well of finite, astrophysically reasonable size. 
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Uacauet &; Krumholzl f|2011h . The instability occurs 
when an interface forms in a gravitational field between a 
low-density radiation pressure-dominated medium on the 
bottom and a higher density, less radiatively-dominated 
medium on top. We can calculate the instability's growth 
rate o f and most unstable w avelength in the linear regime 
using Uacauet fc Krumholzf s formalism. The dispersion 
relation for the instability in the adiabatic, local limit 
(their Equation 77) in our non-dimensionalized units be- 
comes 
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Here the wavenumber k is measured in units of h~ l and 
angular frequency u is measured in units of <?/c s> *. The 
quantities x and /e are the local ratio of radiation pres- 
sure to gas pressure and the local Eddington ratio, re- 
spectively, both computed in the gas immediately above 
the interface. The conditions of adiabaticity and locality 
are satisfied only if (their Equation 82) 
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To evaluate these quantities, we note that before the 
onset of instability F is very close to F everywhere in 
the computational domain, and the values of 6, r , and 
8 g as the base of the layer that goes unstable can be read 
off from the simulation results. Depending on the exact 
time we choose to examine, we find that in run T03F0.05, 
b rs 0.1 — 0.15 and r rj Q g w 1.8. The corresponding fig- 
ures in runs T10F0.25 and T10F0.50 are b « 0.15 - 0.2, 
6 r « Q g w 2.7 and 6 w 0.3, 6 r « 6 S « 2.7. Un- 



fortunately inserting these quantities into Equation ([5t 
indicates that the inequality is not satisfied, indicating 
that we cannot regard the modes as adiabatic, and that 
diffusion is significant. Nonetheless, we can still use 
the dispersion relation (|44[) to obtain an upper limit on 
the growth rate. Inserting these figures into Equation 
(|44[) and numerically solving for w, we find character- 
istic growth rates Im(w) ~ 1, indicating that, if not 
damped by diffusion, the initial seed perturbation we in- 
sert should amplify on a timescale comparable to c S))f /g. 



In practice we see that growth is significantly slower than 
this, almost c ertainly as a result of diffus ive damping. 
Unfortunately Uacauet fc Krumholzl (|2011[) were unable 
to obtain an analytic estimate of the instability growth 
rate in the regime where diffusive damping is significant. 

We can also ask whether the instability we observe in 
our simulations might correspond to any of the other 
types of radiation-hydrodynamic instability that have 
been described in the astrophysical literature. A number 
of authors have studied instabilities in radiation pressure- 
driven flows in the context of clouds near quasars 
(Mat hews fc BlumenthaH Il977t iKrolikl Il977t iMathewsl 
1986). These are quite different than the situation we 
consider in that the opacity comes from resonant ab- 
sorption of ionizing photons, which is therefore linked 
to the recombination rate and hence the density of the 
gas; instabilities arise due to this coupling. Clearly that 
is not the case for the situation we consider, since dust 
opacity is to good appr oximation density- independent. 
iBlaes fc Socrates! (|2003f ) conduct a general analysis of 
local radiative instabilities in optically thick media ap- 
plicable to a wide variety of environments. They find 
that local instability occurs only when there is a mag- 
netic field present or when the opacity contains an ex- 
plicit density dependence, similar to that which applies 
in the quasar case. Our simulations meet neither con- 
dition, and, indeed, t he R RT instability described by 
Uacquet fc Krumho lz (2011) that appears to take place 
in our simulations is an interface instability rather than 
a local one. In the real ISM the dust opacity is, as noted 
before, density-independent. However, the real ISM docs 
contain magnetic fields, and it is conceivable that adding 
a magnetic field to our simulations would produce an 
additional local instability on top of or instead of the in- 
terface one that we find. We discuss this issue further in 
Section 15.5.11 

Perhaps the closest analo g in the l iterat ure to our sit- 
uation is that analyzed by IShaviv] (|2001| ). who studies 
instabilities in atmospheres with a constant, Thomson, 
opacity. He finds that such instabilities can arise when 
the gas is near the Eddington limit, with the instabil- 
ity growth rate dependi ng on the boundary conditions 
(|Blaes fc Socrates! [2003) . The primary difference be- 
tween our work and his is the radiation temperature- 
dependence of the dust opacity of the ISM. This leads 
to an Eddington ratio that is non-constant with height, 
which in turn means that our instability tends to take the 
form of gas at the bottom of the atmosphere where the 
radiation temperature and opacity are high being driven 
upward into a thin shell that collides with gas at larger 
altitudes that experiences a lower radiation temperature 
and opacity. It is this behavior that produces an interface 
and an interface instability. In the absence of radiation 
temperature-dependence in the opacity, radiative accel- 
erations are height- independent, and thin shells should 
not form . In thi s conte xt the global, non-interface insta- 
bility of IShaviv! (|2001l) can occur. The non-linear out- 
come of the two instabilities are likely to be different as 
well. In the case of an ISM with dust opacity, turbulence 
is driven by the circulation of material falling deep into 
the atmosphere, encountering a hot radiation field, find- 
ing itself super-Eddington, and then being blasted up- 
ward; once at height the radiation field pushing on the 
gas is at lower temperatures, the gas falls back, and the 
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cycle repeats. This mechanism obviously cannot operate 
for a gray opacity. 

5.2. Why is Radiative Trapping Ineffective? 

The relatively low values of (/e) and /t rap we find in 
our simulations are surprising. For a laminar matter dis- 
tribution, even if all the gas were at a temperature T* and 
thus had opacity «r,*, we would have /trap = 7* — 1. Con- 
sidering the higher opacities produced by higher temper- 
atures, we might expect /trap ~ 100, as in run T10F0.02, 
and (/e) 3> 1. In the last column of Table [3] we give 
the value of the midplane optical depth tir calculated 
at the end of the simulation using tir « «(T mp )£. Val- 
ues at the start of the calculation are very similar. In 
the three unstable cases, the value of /trap is a factor of 
~ 3 — 6 smaller than tir, with the largest difference oc- 
curring in the most unstable run. How is it possible to 
have such a small value of /trap with respect to the naive 
estimate? To answer this question, it is helpful to exam- 
ine the distribution of radiation flux. We do so in Figure 
using run T10F0.50 at time t/£* = 50 as an example. 
Other time slices and runs when /t ra p is small give sim- 
ilar results, but this time slice, when the instability is 
non-linear but has not yet dissolved into complete tur- 
bulence, provides a particularly clear illustration. As the 
plot shows, the radiation flux is both highly non-uniform 
and strongly anti-correlated with the matter distribution. 
Within the fingers of gas projecting downward, the flux 
is ~ O.li^j while in the narrow channels between the 
fingers the flux is ~ IOFq. This explains how there can 
be so little radiation force exerted on the matter: the 
radiation flux is highly concentrated in low-density, low- 
optical depth channels that contain little mass, so the 
effective optical depth is much less than the true optical 
depth. 

It is important to note, however, that this does not 
mean that the system is effectively optically thin, or that 
an external observer would be able to see directly down 
to the radiation source with little interfering material. As 
the plot shows, even along the vertical paths through the 
simulation domain with the lowest column densities and 
optical depths, corresponding to the channels through 
which most of the flux is focused, the Rosseland mean 
optical depth is always at least ~ 5, and the column 
density is always at least ~ 0.2S. The significant point 
is not that the value of ~ 5 is particularly special; it is 
that the effective momentum imparted by the radiation 
field to the gas can be reduced compared to an estimate 
based on tir without there being optically thin chan- 
nels that would allow direct optical observation of the 
stars providing the radiation flux. The absence of such 
transparent channels does not imply that radiation-gas 
coupling is efficient. 

While the channeling of the radiation is easiest to see 
early in the development of the instability, as shown in 
Figure [9l it continues at later times as well. Figure 
1101 shows the density, temperature, velocity, and radi- 
ation flux distribution in T10F0.50 at the last time slice, 
t = 215£». At this point the radiation pressure-driven 
turbulence is fully developed, and (/e) ~ 1. As the plot 
shows, the radiation flux continues to be highly non- 
uniform. The structure of this atmosphere and its de- 
pendence on model parameters will be the subject of a 
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Fig. 9. — The two images show the gas density distribution (left) 
and magnitude of the radiation flux F (right) at time i/t* = 50 
in run TlOFO.50. The line plots above them show the surface 
density T,(x) and Rosseland mean optical depth tr(x) computed 
along vertical paths at a fixed horizontal position x, i.e. S(x) = 
J„°° p(x, z) dz and tr(x) = J„°° p(x, z)kr(x, z) dz. 




Fig. 10. — Density distribution and Mach number vectors (top) 
and temperature distribution and flux vectors (bottom) in the last 
timeslice (t/t* = 215) of run T10F0.50, showing a snapshot of the 
atmosphere's structure after (/e) — > !• 



future paper. 
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5.3. The Origin of Turbulence in ULIRGs and Dense 

Clusters 

We find that radiation pressure-driven instabilities can 
produce significant turbulence, and thus the hypothesis 
that radiation pressure might be responsible for at least 
some turbulence in dense protoclusters and in the disks of 
ULIRGs seems to be valid. The typi cal Mach numbers i n 
Galactic protoclusters are ~ 10 (e.g. lShirlev et al.l l2003). 
comparable to what we obtain in the simulations, so radi- 
ation pressure-driven instabilities may fully explain the 
turbulence there. (Of course we cannot rule out that 
other effects might contribute as well.) For ULIRGs, 
however, the problem is harder. Even in the most unsta- 
ble run we consider, T10F0.50, the Mach number is still 
only ~ 10. This is roughly an order of magnitude less 
than the Mach numbers seen in real ULIRGs and those 
needed to maintain a Toomre Q parameter near unity 
(jThompson et alJl2005T ). 

What then does this imply about the origin of the 
turbulence in ULIRGs? One possibility is that, given 
the observational uncertainties, ULIRG surface densities 
have been overestimated or ULIRG luminosities underes- 
timated, and in fact they have /e,* ~ 1 over large areas. 
In this case it seems likely that radiation pressure could 
be responsible for driving the turbulence, although the 
problem of how this process picks out Q — 1 would still 
exist. A closely related possibility that answers the ques- 
tion of how Q = 1 is maintained is that the turbulence is 
driven by a limit cycle whereby episodes of star forma- 
tion produce brief periods where /e.* > 1 and the bulk 
of the gas does begin to be ejected. This halts star for- 
mation, and after ~ 4 Myr when massive stars begin to 
die, /e,* drops below unity and the gas falls back. This 
ejection could either be over a large area of the disk, 
or it coul d occur locally in wi thin forming clusters. For 
example, iMurrav et all ()2010D suggest a cycle in which 
Toomre-mass clumps repeatedly form, undergo gravita- 
tional collapse, and then are disrupted by star clusters 
that rapidly reach luminosities such that /e,* > 1. The 
ejecta from these disrupting clusters provide the turbu- 
lence. 

Alternately, the turbulence in ULIRGs might not 
be radiatively driven at all. Instead, it might be 
driven by gravitational instab ilities in ULIR G disks 
(jKrumholz & BurkertlfeOlOl : iForbes et alJl2012t ). and be 
maintained by the energy of inward gas mass through the 
galactic potential rather than by stellar feedback. The 
turbulence could als o be driven by c o smic ray rather than 
radiation pressure (Ipavic hl 119751: [B reitsch werdt et all 
ll991ilJubelgas et alJl2008l : ISocrates et alJl200^ T 

5.4. Implications for Analytic and Numerical Models of 
Radiation Pressure Feedback 

Our work suggests that radiation pressure alone is un- 
likely to be able to eject much matter or drive signif- 
icant winds when K,(T mp )pF/c < Eg < K,(T c d gc ) pF / c, 
i.e. the regime where the radiation flux is sub-Eddington 
at the top of the disk where the temperature is low, 
but super-Eddington at the midplane where the tem- 
perature is high. In this regime, we find that even a 
very laminar initial density field will spontaneously re- 
arrange itself, via radiation Rayleigh- Taylor instability, 
into a configuration whereby the mean Eddington ra- 



tio is unity or slightly smaller. The final-state struc- 
ture we obtain for T10F0.5 is shown in Figure [10J In 
this limit a small amount of mass could conceivably be 
blown off as a wind, though this does not occur in our 
simulations, but it is not possible to eject enough mat- 
ter to materially reduce the star formation efficiency. 
Bulk ejection appears likely only when the radiation 
field exerts enough force to be super-Eddington even 
using the lower opacity present at the top of the at- 
mosphere. This condition is expected to be met un- 
der many circumstances, particularly in young star clus- 
ters ([T homps on et all 120051: iKrumholz fc Matznerl 120091 : 
iFall et al.ll2010t IMurrav et al.ll2010fl . and thus radiation 
pressure may still be a significant factor in limiting star 
formation efficiencies. However, models that rely on a 
large enhancement of the radiation force due to higher 
opaci ties in the warmer regions (e.g. IMurrav et al.ll2010l 
12011ft may need to be reevaluated. In particular, in 
the three unstable models considered here, the ratio of 
tir ~ K(T mp )S to / trap is ^ 3 - 6 (Tabl e |3). Sim- 
ilarly, the mode ls of IKrumholz fc Matznerl (|2009f) and 
IFall et all (|2010t l may also need to be reconsidered since 
they assume /trap ~ 1- We find that /trap can signif- 
icantly exceed unity, which suggests that, if radiation 
does overcome gravity and eject matter, these models 
might underestimate the total momentum input to the 
gas by a factor as large as ~ 5 — 40. However, since none 
of our runs successfully launch a wind, the question of 
how much momentum will be transferred to the gas if 
(/e) > 1 and a wind is launched is not yet settled. 

For subgrid models of radiative pressure feedback in 
numerical simulations, the main implication of our work 
is that the true value of /trap is likely to be somewhat 
lower than tir, and that /trap will not increase in direct 
proportion to the total gas column density, or some proxy 
for it. Instead, we find that a reasonable estimate is that 
when /e,* > /E.crit, the system will adjust to give 
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Note that with this scaling /t rap still increases in direct 
proportion to E, as assumed in many models, but that 
there is an additional dependence on the flux, the opacity, 
and the gravitational force holding the gas together that 
ensures that the mean Eddington ratio (/e) saturates 
at unity. Simulations based on subgrid models in which 
/trap is allowed to reach values such that (/e) is larger 
than unity need to be recomputed with lower values of 
/trap to check which results are robust. 

It is illuminating to compare our resu lt to both 
the fidu cial estimate of /trap adopted by Hopkins et al. 
(2011g), and to the alternative models presented in Ap- 
pendix B of their paper, where th ey attempt to take in to 
account photon leakage (see also IMurrav et al.l J2010D1. 
As noted above, in their fiducial models Hopkins et al. 
take /trap = tir — 1- In the alternative models, they 
consider media with a variety of column density distri- 
butions motivated from observations of molecular clouds 
and simulations of turbulence. For the case that pro- 
duces the most deviation from their fiducial estimate, 
that of an exponential distribution of column densities 
with a powerlaw tail, they find that for large optical 
depth the effective value of /trap approaches /trap ~ 
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^/(r/2r[l/cr]r II ( — 1, where T is the T function, ttr 
here is the optical depth that the medium would have 
if it were uniform, and a is th e standard dev iation of 
the column density distribution. Hopkins et al.l consider 
values of a from 0.2 — 1.0, and for the largest value of 
a they consider, at high tir this prescription therefore 
reduces to ,f t ra P ~ >At r/2 - 1. Thus, for / trap > 1, the 
ratio of iHopkins et all s estimates of /t rap to the upper 
limit we measure in our simulations (Equation I51[) is 
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where the first quantity is for Hopkins et all s fiducial 
estimate, and the second is for their alternative model 
with cr=l. Similar estimates for / t rap in a turbulent 
medi um were ma de by iMurrav et all l|2010l )). 

Hopkins et all s fiducial estimate for tir is Tir = 
5(E/g cm -2 ), whereas, for a normal galaxy in which 
T* ~ 40 K, we will have r* ~ Ttr/10. Taking this as 
a rough estimate, Equation (1521 becomes 



/1 



trap, Hopkins 



/( 



trap,lim 



10/) 



E,* or 



(53) 



Thus we find that the IHopkins et al.l (|2011bT ) estimate 
of /trap exceeds the limit we measure whenever /e,* > 
0.1 for their fiducial estimate, or when /e,* > \/50/tir 
for their lowest alternative estimate. The first condition 
is likely met in many places in their simulations. The 
second is somewhat harder to satisfy, but it still likely 
to be met at the highest optical depth locations in their 
simulations. 

The prim ary reason IHopkins et al.l and 
IMurrav et al.l estimate larger effective trapping fac- 
tors than we measure is that while their models consider 
the possibility that the gas might be non-uniform, 
they do not consider that the radiation might also 
be non-uniform, and that its non-uniformity might 
be correlated with that of the gas. Figures [9] and ITOl 
show that, once the instability is fully established, this 
radiation-matter correlation is an essential feature of 
the radiation-gas coupling that these analytic models do 
not capture. 

5.5. Caveats and Future Work 
5.5.1. Magnetic Fields 

It is important to point out some potential limi- 
tations of our results, which suggest avenues for fu- 
ture investigation. The first is that we have omitted 
magnetic fields, which are clearly present in the real 
interstellar medium. Simulations of mechanical feed- 
back f rom both protq s tellar o utflows dLi fc Nakamural 
[20M INakamura fcH 120071 : IWang et al.l 12010ft and 
photoionizati on-driven "champagne" flows from molec- 
ular clouds (|Gendelev k, Krumhold I2012D show that 
strong, ordered magnetic fields can significantly enhance 
the effects of feedback by providing a mechanism to 
transfer momentum between parcels of gas, and thus dis- 
tribute energy and momentum more broadly through a 
gas cloud. This effect could conceivably operate here, 
and raise (/ E ) and / trap . 

However, the effect is only likely to be important 
if the magnetic fields are dynamically strong; other- 



wise field lines will bend rather than exert significant 
forces. In ULIRGs, indirect arguments based on obser- 
vations of synchotron emission suggest that the mag- 
netic energy density is sub-equipartit ion compared to 
the turbulent or gravitational energies ([Thompson et all 
l2009HLacki et alJl2010HBateiat et al.ll2011f l. which sug- 
gests that magnetic effects are unlikely to be important 
there. The situation is less certain for massive star clus- 
ters in normal galaxies like the Milky Way. In star- 
forming clouds in the Milky Way and nearby galaxies, 
Zeeman and polarization observations indicate that there 
are ordered magnetic fields i n rough equipartition with 
turbulence dCrutcherl fl999t iTroland fc Crutcheii [20081: 
Chapman et al. 20111 : iLi fc H cnning 20 111 however see 



Padoan et al.ll2004l for a contrasting view) , in which case 



magnetic effects could conceivably alter (/e) and /trap- 
However, all of the regions studied thus far are far from 
the regime of massive clusters forming at high volume 
and column density where radiative forces might be im- 
portant, and it is unknown if the magnetic fields in such 
regions are also in equipartition. In any event, in future 
work it would be useful to investigate whether the inclu- 
sion of magnetic fields allows larger values of (/e) and 
/trap, and, if so, how large the field must be to achieve 
this effect. 

Another possible effect of magnetic fields is that they 
might make the gas subject to the loca l photon bub- 
ble instability of iBlaes fc Socrates! (|2003l ). as discussed 
in Section [5. II In this case a local instability might exist 
on top of, or in place of, the interface one that occurs in 
our non-magnetic simulations. It seems unlikely that this 
additional instability would strengthen matter-radiation 
coupling, but it is conceivable that it could weaken it 
even further relative to what we have found. 

5.5.2. Planar Versus Spherical Geometry 

For simplicity we have chosen to examine a planar ge- 
ometry, since a spherical geometry would introduce a 
third parameter (which we could take to be the radius 
of curvature measured in units of /i*) to our description 
of the system. The main disadvantage of our planar ap- 
proach is that in planar geometry the escape speed is not 
well-defined, and this precludes one possible mechanism 
for launching a wind. We find that, in steady state, (/e) 
reaches unity. However, there is a transient phase before 
the instability develops when (/e) is larger than unity. In 
planar geometry this makes little difference, because the 
momentum the gas absorbs during this transient is in- 
sufficient to escape, and the matter will always fall back. 
For a spherical geometry, however, it is conceivable that 
the gas could reach escape speed during the initial tran- 
sient when (/e) > 1, allowing it to escape even though 
(/e) falls back to unity once steady state is reached. 

This is unlikely to be a significant effect for the pa- 
rameter regime we have explored. For galactic disks, 
non-spherical effects become significant only on scales 
comparable to the radial scale length, which is of or- 
der kpc. For comparison, for the ULIRG-like scalings 
shown in Table [1] all our simulations boxes arc < 1 pc 
in size. Thus we can be confident that non-planar effects 
are unimportant in the galactic case unless the Edding- 
ton ratio, and thus the typical height to which matter 
is carried by radiation, is far larger than in the parame- 
ter regime we've explored. For the case of individual star 
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clusters, which can be ~ 1 pc in size, the highest Edding- 
ton ratio case may be marginally in the regime where 
curvature effects become significant. However, even in 
this case most matter only reaches ~ 0.3 pc before turn- 
ing around. Moreover, in a real star cluster where the 
initial matter distribution is not laminar, and where the 
radiation flux rises smoothly rather than turning on in- 
stantly as in our simulations, the amount of momentum 
deposited during the initial transient is likely to be less 
than in our simulations. 

Thus we conclude that curvature effects are unimpor- 
tant for the parameter regime we have explored, and are 
only likely to become important for systems with Ed- 
dington ratios significantly above 0.5. 

5.5.3. 2D versus 3D 

For reasons of computational efficiency we have 
conducted two-dimensional rather than fully three- 
dimensional simulations. While the dimensionality does 
not change the gr owth rate of the RRT inst ability in 
the linear regime (jJacquet fc Krumholzl l201lT ). it may 
affect the non-linear growth rate and fully saturated 
state of the instability. Even for simple fluid Rayleigh- 
Taylor instability the relationship between two- and 
three-dimensional results is still not fully understood, 
though numerical results indicate the non-linear growth 
rate is a factor of ~ 2 faster in the 3D case (e.g. 
I Young et al.l 1200 ll and references therein). No compa- 
rable study exists for RRT, nor have we conducted such 
numerical experiments. However, it seems unlikely that 
the difference between 2D and 3D will be astrophysically 
important. The time required for the instability to reach 
full non-linear saturation in our simulations is smaller 
than essentially any other timescale relevant on galactic 
or star cluster scales, and a factor of order unity change 
in the non-linear growth rate would not alter this. More- 
over, the non-linear saturated state we obtain is defined, 
as is the the case for many saturated instabilities, by the 
system self-regulating to a state of marginal stability, 
(/e) ~ 1. It also seems unlikely that this self-regulation 
will fail in 3D. Thus the differences between 2D and 3D 
are unlikely to be astrophysically significant. 

5.5.4. Dependence on Choice of Initial Conditions and 
Simulation Box Size 

Our computations use an idealized initial setup and 
computational box size chosen to allow us to follow the 
growth of the instability through its linear and then non- 
linear development, while ensuring that the smallest im- 
portant size scale /i* was always at least marginally re- 
solved. It is fair to ask how the results might change 
for a more realistic setup that would correspond more 
closely to a real star cluster or ULIRG. In considering 
this question, it is helpful to separate the question of the 
initial horizontal and vertical structures. 

In the horizontal direction, our initial conditions are 
characterized by a fairly small amplitude initial pertur- 
bation with power only at a single horizontal wavelength 
whose physical size is very small compared to real astro- 
physical systems - for example, run T10F0.25 has only a 
25% density perturbation at a physical size of 0.082 pc. 
Moreover, in the horizontal direction our entire computa- 
tional box is only a factor of 2 larger than this. Since real 



protoclusters and ULIRGs are both supersonically turbu- 
lent and much larger than this, a real system would have 
much larger amplitude perturbations on much larger hor- 
izontal size scales. How would this affect our results? 

In the linear regime, the analytic treatment discussed 
in Section 15.11 is a helpful guide. The linear analysis 
shows that modes with short horizontal wavelengths are 
damped by horizontal diffusion of the radiation, which 
reduces the linear growth rate of the instability. Pertur- 
bations with longer horizontal wavelengths would be less 
damped, and thus should grow faster, with the growth 
timescale approaching t* ~ 1 at long wavelengths, rather 
than £» ~ 10 — 100 as found in our simulations^ More- 
over, with larger initial perturbations, fewer e-foldings 
would be required to reach the non-linear regime. Thus 
we conclude that in a real astrophysical system, the tran- 
sition to the non-linear regime is likely to be significantly 
faster than in our simulations. Once in the non-linear 
regime, it seems likely that the same argument we made 
in considering 2D versus 3D applies: the non-linear state 
is characterized by the system self-regulating to a state 
of marginal stability, and it seems unlikely that this self- 
regulation depends on the size of the computational box 
or the initial state. Thus in the saturated state we do not 
expect the result that (/e) ~ 1 would be different. How- 
ever, the density structure might well extend to larger 
size scales in a larger computational domain. For ordi- 
nary fluid RT instability, in the non-linear regime the 
dominant mode is always on the largest scale permitted 
by the computational domain or experimental apparatus, 
and since RRT instability begins to behave like ordinary 
RT instability at very long wavelengths, this would likely 
be true for our problem too. 

Our initial state is also very thin in the vertical direc- 
tion, since it corresponds to a gas layer supported only 
by thermal pressure at constant temperature. In a real 
ULIRG or protocluster, the turbulent velocity greatly ex- 
ceeds the thermal sound speed, and puffs up the gas to a 
much larger height, which is comparable to or larger than 
the vertical height we reach at the end of our simulations 
when the RRT instability is fully saturated. If we started 
with such a turbulent, vertically extended state but did 
not include any mechanism to drive the turbulence other 
than the RRT instability itself, then it seems likely that 
the result would be that the turbulence would decay un- 
til the gas reached the velocity dispersion dictated by 
the instability. The timescale for this to occur would 
likely simply be the turbulent decay time, which would 
be comparable to the crossing time in the initial verti- 
cal distribution. On the other hand, if there were some 
driving to maintain the turbulence, it is likely that this 
would dictate the gas structure. This probably could not 
cause (/e) to exceed unity, since this would imply that 
radiation rather than the other mechanism had become 
dominant. However, extra structure in the gas produced 
by turbulent driving could cause gas-radiation coupling 
to change in unexpected ways, and might conceivably 
lead to either stronger or weaker coupling of radiation 
and matter than would have occurred without the extra 

8 Conversely, shorter wavelength perturbations than those we 
have used would grow more slowly or not at all; indeed, in some low 
resolution simulations we conducted using even smaller horizontal 
box sizes, we did not see development of RRT instability. However, 
no real astrophysical system is this small horizontally. 
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driving. 



5.5.5. Nature of the Radiation Source 



Another caveat is that we have explored the effects of a 
radiation source that is steady, uniform, and located just 
below a gas layer, whereas in reality the young stars that 
provide radiation pressure in a young cluster or a ULIRG 
are spatially clustered, time- variable (as new stars form 
and old ones die), and are mixed with the gas. It seems 
highly unlikely that any of these effects would increase 
(/e) or /trap, since making the radiation field even less 
uniform would only further weaken matter-radiation cou- 
pling. On the other hand, as noted in Section 15. 3| it 
seems possible that such non-uniform radiation source 
could produce larger turbulent velocities than we find in 
our simulations with a uniform radiation source. Again, 
followup simulations are needed to check this effect. 

5.5.6. Numerical Treatment of the Radiative Transfer 
Problem 

A final concern is the quality of the two-temperature 
flux-limited diffusion (FLD) approximation we use to 
treat the radiative transfer problem. One significant 
concern is our treatment of the direct (as opposed to 
dust-reprocessed) radiation field. The direct radiation 
field possesses two properties that are not well-captured 
by the FLD approximation. First, its spectrum is at 
much higher frequencies than would be predicted by our 
blackbody approximation; this is significant because dust 
opacities increase strongly with frequency. Second, the 
radiation field is highly directional, as opposed to the dif- 
fuse field assumed in the FLD approximation. In the case 
of the formation of a sing le massive star from an initially 
laminar protostellar core, iKuiper et al.1 (|2011h show that 
omitting the direct radiation field can lead to an underes- 
timate of the expansion rate of radiation pressure-driven 
cavities, and that this in turn can affect whether a cavity 
goes RRT unstable before it blows out of its parent core. 

This effect is unlikely to be important for the simu- 
lations we report here, simply because, by construction, 
the direct radiation pressure force is unimportant. The 
strength of the direct radiation pressure force relative 
to gravity can be characterized by the mean Eddington 
ratio considering only direct photons, which is given by 
Equation (J43I) evaluated with /trap = 0: 



If X /e,* 

UE.dir/ — ■ 



(54) 



This gives (/ E dir) = 0.002,0.17,0.025, 0.05, and 0.09 
for runs T10F0.02, T03F0.50, T10F0.25, T10F0.50, and 
T10F0.90 respectively. Thus we see that the direct ra- 
diation force is an order of magnitude or more weaker 
than gravity in our runs, and its omission is therefore 
not likely to produce significant effects. 

It is important to note that the relative unimpor- 
tance of the direct radiation field compared to gravity 
in our simulations is a direct result of our parameter 
choices, which are somewhat different from those appro- 
priate to the s ingle massive protostellar cores studied by 
IKuiper et al.l (|2011l) . To see why these systems are in a 
somewhat different regime than ULIRGs or star clusters, 
it is helpful to re-express (/E.dir) in terms of dimensional 
parameters: (/E,dir) = Fo/Y*gc = L/Mgc, where in the 



second step we have multiplied through by a fiducial area 
to turn Fq into a luminosity L and £ into a mass M. 
If the gravitational field comes predominantly from the 
self-gravity of the object question, then g « GE, and we 
have 



(/E,dir) ~/i 



L/M* 
GT,c 



--4.8 x 10" 



1/2 



(L/M: 



*Jo^o 



(55) 



where (L/M,)o = (L/RQ/(L Q /M Q ), S = E/l g cm" 2 , 
and for convenience we have defined /* as the stellar mass 
fraction, i.e. the fraction of the object's mass that is in 
luminous stars rather than gas. Note that this equation 
for the importance of direct radiation forces is, up to 
factors of order unity, the same as the result derived by 
iFall et all (|2010l their Equation 5). Individual massive 
protostellar cores, massive protoclusters, giant clumps in 
high redshift galaxies, and ULIRG disks all have similar 
surface densities £ sw 1 g cm" 2 , but they have very dif- 
ferent stellar light to mass ratios. For example, a 100 
Mq star on the zero-age main seque nce ( ZAMS) has 
L/M„ = 1.3 x 1O 4 L O /M (|Tout et al.M1996h . so (/ E ,dir) 
can easily be significantly greater than unity, particularly 
in a rotating system where rotational flattening lowers 
the surface density in certain directions. In contrast, 
a cluster of ZAMS stars drawn fr om a fully sampled 
IMF will have L/M* « 1O 3 L Q /M (jMurrav fc Rahman! 
I201CO . an order of magnitude lower. Since L/M* declines 
as a stellar population ages, (/E,dir) will likely be smaller 
than unity in ULIRG disks or in the giant clumps in high 
redsh ift galaxies with dynamic al times longer than ~ 4 
Myr (|Krumholz fc Dekell l2010l ) . although (/ E ,dir) might 
still exceed unity in subregions where the stellar popula- 
tion is predominantly young. We can therefore conclude 
that, even if the direct radiation field can be dominant 
for single massive stars and cores, it is at most an order 
uni ty effect for massive star clusters and ULIRGs (but 
see iMurrav et all (J2010) and iKrumholz fc Dekell (J2010T ) 
for more discussion) 



Furthermore, as iSocrates et al.1 (J2008I ) point out, the 
above calculation may overestimate the importance of 
the direct radiation force if the geometry of the emitters 
is uniform, because it does not consider the fact that the 
direct radiation forces supplied by stars can cancel. A 
single massive star supplies a radially outward radiation 
field, but the many stars in a cluster or disk will push the 
gas with which they are intermingled in many different 
directions, and thus there will be some cancellation, with 
the exact amount d epending on the geometry. However, 
note that Hopkins et al. (2011b) argue that the effects 
of cancellation are not strong in simulations of radiation 
pressure feedback in simulations of whole galaxies. 

Even if omission of the direct radiation field is not 
problematic in our physical situation, one can still worry 
about other aspects of the FLD approximation. FLD 
requires that the region of interest be optically thick, 
and our simulation domain certainly satisfies that crite- 
rion: tir exceeds 10 in all the runs and exceeds 100 in 
three of them (see Table [3]). However, since FLD dis- 
cards information on the directionality of the radiation 
field, it underestimates beaming of the radiation field in 
localized patches that may be optically thin. Quantify- 
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ing the net effect of this error is difficult, and generally 
requires direct comparison of FLD results with more ac- 
curate methods such as variable Eddington tensor, dis- 
crete ordinates, or Monte Carlo. One such comparison, 
done in the context of a r adiation-dominated accretion 
disk bv iDavis et al.1 (J2012J ). suggests that FLD tends to 
overestimate the vertical force exerted by the radiation, 
because it underestimates the ability of the radiation to 
stream in the horizontal direction. If this were true of 
our problem it would suggest th at we have ov erestimated 
/trap- However, we caution that IDavis et al.l 's problem is 
not completely an alogous to o urs (for example the dom- 
inant opacity for IDavis et al.l is the scattering opacity 
of free electrons, rather than the absorption opacity of 
dust), and it is not clear if the same would be true for 
our case. 

6. CONCLUSIONS 

We present numerical simulations of a strong radiation 
flux passing through a column of gas confined by grav- 
ity. This configuration is a reasonable proxy for a galac- 
tic disk or a section of a young star cluster in which the 
radiation from young, massive stars passes through the 
dusty interstellar medium and exerts dynamically sig- 
nificant forces. This system is characterized by two di- 
mensionless numbers, t» and /e,*- The former describes 
the optical depth of the gas column computed using the 
temperature at its surface, and the latter describes the 
ratio of radiation pressure force to gravitational force in 
this gas. We use these simulations to study whether ra- 
diation is able to drive turbulence or produce winds in 
the regime where the radiation force is sub-Eddington 
for the cool gas at the top of the disk / the edge of the 
young cluster (/e,* < 1), but the gas is optically thick 
(t* > 1), so that the higher temperatures make the gas 
super-Eddington near the disk midplane / cluster cen- 
ter. The disks of ULIRGs and some young star clusters 
in non-ULIRG galaxies are in this regime. 

We find that, in this regime, radiation forces drive the 
matter into a thin sheet which then breaks up due to ra- 
diation Rayleigh- Taylor instability (Figs. EJ El [TO]). Once 
this instability reaches its full non-linear development, 
it taps a fraction of the radiation energy and uses it to 
drive supersonic turbulent motions in the gas. In the pa- 
rameter regime appropriate to ULIRGs and young clus- 
ter clusters the velocity dispersion can approach Mach 
10, sufficient to fully explain the turbulence observed in 
young protocluster gas clouds in the Milky Way. ULIRGs 
show significantly greater velocity dispersions, which sug- 
gests that either radiation pressure-driven instabilities 



cannot drive turbulence at the required levels, or that 
ULIRGs are closer to the Eddington limit at the cool 
tops of their atmospheres than current observations sug- 
gest; this is entirely possible, given the observational un- 
certainties. 

We also find that radiation Rayleigh- Taylor instabil- 
ity leads to a configuration in which the matter is con- 
centrated in filaments, while the radiation flux is con- 
centrated in low-density channels (Figs. HI [TO]) . In this 
configuration radiation is not fully trapped by the gas. 
As a result, the actual mean force applied to the gas 
never exceeds that applied by gravity, despite the fact 
that gas near the midplane is super-Eddington. We find 
radiation passing through slabs of matter with /e.* < 1 




Fig. 11. — Comparison of the total velocity dispersion a (top 
panel), mean Eddington ratio (/e) (bottom panel, left axis), and 
trapping factor /trap (right axis) in runs T10F0.25 (thick lines) and 
T10F0.25-LR (thin lines). The results are clearly quite similar. 
To avoid clutter, we do not show a x , <r z , or the trapping factor 
considering only material with v z > 0, as we do in Figures [7] and 
[8] but these lines are also very similar in the two runs. 

does not drive significant mass ejection or a noticeable 
wind even if t*/e,* S> 1. The ratio of the momentum 
transferred to the gas to that carried by the radiation 
field, defined as 1 + /trap where /trap is the trapping fac- 
tor, reaches a value (given by Equation I5ip that ensures 
that the mean Eddington ratio is unity. These numerical 
results are in conflict with the assumptions built into an- 
alytic models and sub-grid implementations of radiation 
pressure feedback in numerical simulations, which either 
limit /trap ~ 1 or allow substantially larger values of 
/trap- Simulations based on these assumptions should be 
recomputed using our improved determination of /trap- 



APPENDIX 

A. RESOLUTION STUDY 

To check the dependence of our results on the numerical resolution of the simulations, we rerun simulation T10F0.25 
at half resolution; we call this run T10F0.25_LR, and describe its properties in Tabled) Since the instability in its 
fully saturated state is chaotic, we do not expect the results of runs T10F0.25 and T10F0.25_LR to be identical in 
more than a statistical sense in the non-linear regime. Figure [TT1 shows a comparison of the time evolution of the gas 
velocity dispersion, Eddington ratio, and trapping factor in the two runs, and Table [3] gives the quantitative results. 
As the Figures and Table show, the results are in line with what one would expect for a converged calculation. At early 
times, before the instability becomes non-linear, er, (/e), and /trap are nearly identical in the two runs. As time goes 
on, the instability goes non-linear, and the flow becomes chaotic, the two runs diverge, but they remain statistically 
nearly identical. In particular, the primary result that the mean Eddington ratio self-regulates to unity at late times 
once the instability is fully saturated is found at both resolutions. The time-averaged values of a and ttr, differ by 
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~ 10% between the two runs, but these differences are smaller than the fluctuations in time in these quantities found 
in each run, and thus are consistent with simply being the results of random sampling of the chaotic flow pattern. 
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